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Boussinesq approximation of the Cahn-Hilliard-Navier-Stokes equations 



O 

(N 

c/3 



V 

i 

o 

>v 
^ : 

Oh. 



> 

in 

oo 
q 

ON 

o 
o 



X 



Anatoliy Vorobe\EI 
Energy Technology Research Group, 
School of Engineering Sciences, 
University of Southampton, SO 17 IB J, UK 
(Dated: September 7, 2010) 

We study the interactions between the thermodynamic transition and hydrodynamic flows which 
would characterise a thermo- and hydro-dynamic evolution of a binary mixture in a dissolu- 
tion/nucleation process. The primary attention is given to the slow dissolution dynamics. The 
Cahn-Hilliard approach is used to model the behaviour of evolving and diffusing interfaces. An im- 
portant peculiarity of the full Cahn-Hilliard-Navier-Stokes equations is the use of the full continuity 
equation required even for a binary mixture of incompressible liquids, firstly, due to dependence of 
mixture density on concentration and, secondly, due to strong concentration gradients at liquids' 
interfaces. Using the multiple-scale method we separate the physical processes occurring on different 
time scales and, ultimately, provide a strict derivation of the Boussinesq approximation for the Cahn- 
Hilliard-Navier-Stokes equations. This approximation forms a universal theoretical model that can 
be further employed for a thermo/hydro-dynamic analysis of the multiphase systems with strongly 
evolving and diffusing interfacial boundaries, i.e. for the processes involving dissolution/nucleation, 
evaporation/condensation, solidification/melting, polymerisation, etc. 

PACS numbers: 47.51. +a, 47.61. Jd, 68.05.-n 



I. INTRODUCTION 

Fluid mixtures can be classified as being immiscible 
(e.g. oil/water mixture) or miscible (e.g. honey/water 
mixture or mixture of two gases). An immiscible system 
is characterised by a strict interfacial boundary. This 
boundary cannot be crossed by the molecules of adjoining 
liquids, which, on molecular scale, can be explained by 
a high potential barrier due to different intermolecular 
interactions within the mixture components. At macro- 
scale, the coefficient of surface tension is introduced to 
define the macroscopic effects of the molecular potential 
barrier. 

The mixture of two gases is an example of a directly 
opposite case. As intermolecular forces between gas 
molecules are negligibly small, no potential barrier at the 
gases' boundary exists, and gas molecules of initially sep- 
arated components freely co-diffuse. So there is no sense 
in introduction of an interface, and hence, there is no 
surface tension on the gases' boundary. 

The focus of the current work is on the miscible mix- 
tures of two liquids, for which the intermolecular forces 
cannot be neglected. As an everyday example, the be- 
haviour of a droplet of honey in water may be considered: 
for such a droplet, a strict interface is visible for a long 
period after immersion of a droplet into water; but after 
a slow dissolution process honey /water mixture becomes 
homogeneous. Again, the existence of a strict interface 
between mixture components, on molecular scale, should 
be associated with a potential barrier. But, in compari- 
son with the immiscible case, for miscible interfaces, some 
molecules with sufficiently high kinetic energies are able 
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to cross this barrier. It may be further assumed that 
the molecular-flux through an interface gradually dimin- 
ishes the barrier's height which would result in growing 
numbers of molecules being diffused from one phase to 
another, and, ultimately, in a complete dissolution of a 
droplet. 

Thus, the concept of interface is required to describe 
the behaviour of a slowly miscible system. The inter- 
face, at macro-scale, will be characterised by the surface 
tension coefficient. In contrast to immiscible systems, 
surface tension coefficient varies as the dissolution pro- 
cess progresses, see e.g. Ref. [l]. The cases of completely 
and partially miscible liquids (e.g. honey/water and bu- 
tanol/water mixtures) would differ from each other by 
the fact that the surface tension coefficient decreases to 
zero in the first case (as the interface disappears) and is 
dynamically variable over dissolution process but always 
remains non-zero in the second case. 

Concluding this discussion, we would like to state that 
the mixing of two (even completely miscible) liquids can- 
not be reduced to a simple diffusion (as it can be done for 
mixing of two gases). The surface tension effects have a 
two- fold influence on the dissolution dynamics, both the 
morphology of the interface and the rate of mass transfer 
through the interface are affected. We may argue that 
the surface tension is sufficiently high to exclude the co- 
diffusion of molecules between immiscible liquids. In the 
case of miscible systems, the mass transfer through the 
interfacial boundary is not zero but its rate is restricted 
by the surface tension effects. 

A first review devoted to the physics of slowly miscible 
systems can be found in Ref. [D. The phase-field ap- 
proach was adopted to model a hydrodynamic evolution 
of a miscible interface. The authors also pointed out an- 
other physical effect that should characterise evolution of 
binary mixtures, namely, the quasi-compressibility of the 
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hydrodynamic equations: even for a mixture of two in- 
compressible liquids, the fluid velocity is a non-solenoidal 
field due to dependence of mixture density on concentra- 
tion. 

The first comprehensive numerical studies for misci- 
ble interfaces were however based on the governing equa- 
tions in which the transport of a second component was 
modelled as an impurity [2|, i.e. with no account for 
the surface tension effects. The Cahn-Hilliard-Navier- 
Stokes equations, that include all essential physical ef- 
fects previously discussed by Joseph and Renardy, have 
been derived in Ref. H. These equations fully define 
the hydrodynamic behaviour of miscible binary mixtures. 
Nevertheless, the recent numerical studies of the misci- 
ble multiphase systems are still based on the impurity- 
like equations, in which however the Korteweg stress has 
been added into the momentum balance . We should 
note that the Korteweg stress is only one of the surface 
tension effects, since the rate of phase change is also de- 
termined by the surface tension effects, which however 
is not taken into account in any recent hydrodynamic 
model for a miscible multiphase system. 

Hydrodynamics of immiscible systems has been stud- 
ied either for density-matched fluids 0-E1, for such flu- 
ids, the quasi-compressibility effects are eliminated, or by 
using the Boussinesq-like approximation [IcHl4j |. Since 
within the phase-field approach it is necessary to assume 
that the density gradients are large at least in some parts 
of a computational domain, this makes justification of 
the Boussinesq approximation difficult. The role of the 
quasi-compressibility effects for a particular system was 
examined in Refs. [l5| (for an analysis of one-dimensional 
diffusion within a pipe) and (for an analysis of mis- 
cible displacements in capillary tubes). Following Ref. 
[HI the velocity field was divided into incompressible and 
expansion parts. In both works, the estimations showed 
that the expansion part of the velocity field is negligi- 
bly small. Nevertheless, even such a statement is not 
sufficient to prove the use of the Boussinesq approxima- 
tion. At present, we are unaware of any paper where 
the Boussinesq approximation of the full Cahn-Hilliard- 
Navier-Stokes equations has been strictly derived. 

To conclude this introductory section, we would like 
to mention that applications that involve the miscible in- 
terfaces are ubiquitous. They include solvent extraction, 
cleaning, removal of oil spills, waste treatment, enhanced 
oil recovery, drug delivery, etc. Frequently, it is assumed 
that the rate of dissolution and the rate of flow change 
are not comparable, e.g. dissolution is a slow process and 
the changes, caused by hydrodynamic flows, happen at 
much faster rate, or, on the contrary, nucleation is a very 
fast process and the hydrodynamic flows are much slower. 
For such cases, different simplified models were proposed 
which allow these processes to be considered separately. 
We however are interested in slow flows when the typical 
dissolution and convective time-scales are comparable; 
an example can serve the hydrodynamic flows in capil- 
lary tubes or in porous media. Other important examples 



where the thermodynamic and hydrodynamic processes 
interact are the flows in near-critical systems. 



II. CAHN-HILLIARD-NAVIER-STOKES- 
EQUATIONS 

A. Phase-field approach 

There are several different approaches used to describe 
the behaviour of a multiphase system. An immiscible in- 
terface has usually very small thickness (of several molec- 
ular layers), which makes the use of the classical Laplace 
approach to model the interface as a surface of discon- 
tinuity (of zero thickness) being well justified. Follow- 
ing this approach, two systems of Navier-Stokes equa- 
tions are separately solved in each phase and so ob- 
tained solutions are matched using boundary conditions. 
This approach is however difficult to apply when interfa- 
cial boundaries undergo complex (topological) modifica- 
tions. For similar problems, the phase-field (also called 
the diffuse-interface) approach becomes more suitable. 

The main idea of the phase-field approach is to arti- 
ficially smear the interfacial boundary and one system 
of Navier-Stokes equations is then solved for the entire 
multiphase system. All variables experience strong but 
continuous changes through the interface. The position 
and shape of the interface is tracked either by introduc- 
tion of additional scalar function (e.g. level-set function 
in the level-set method) or by using a natural variable 
(like concentration field) . A new volume force is also in- 
troduced to mimic the surface tension effects. The fluid 
behaviour in the limit of zero surface thickness is usually 
analysed. We however should notice that the interface 
smearing is physically justified for the systems near ther- 
modynamic critical point. 

In the current work, the phase-field approach is utilised 
to define the evolution of a multiphase system. The mix- 
ture components are called solvent and solute. To char- 
acterise a state of the binary mixture the concentration 
field is used, which we define as the mass fraction of the 
solute in the solvent phase. 

B. Thermodynamic model 

To define a thermodynamic state of a binary mixture 
we need to introduce the free energy function. This func- 
tion can be either written based on experimental data or 
derived from a molecular level theory. Only for illus- 
tration purposes, let us show how the free energy of a 
heterogeneous system should be defined if the classical 
Laplace approach is used. In this approach the phases 
are treated separately and the free energy function in- 
cludes three terms, 

F = V 1 f 1 (p ai ) + V 2 f 2 (p 2) + S(V 1 ), (1) 
S(Vx) = a^Gir) 1 ^^ 3 . (2) 
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FIG. 1. The phase diagram for a binary system with an upper 
critical solution temperature. 



Here, the first and second terms are, respectively, the free 
energies of the first and second components of a mixture. 
The last term, S(Vi), is an interfacial energy traditionally 
defined through the surface tension coefficient a. V\ and 
V% are the volumes of the first and second phases. 

In the phase-field approach, one system of equations 
is used to define the behaviour of an entire multiphase 
system. It is convenient to introduce the specific free 
energy function /. To take the surface tension effects 
into account, Cahn and Hilliard proposed to define 
/ not only as a function of density and concentration but 
also as a function of concentration gradient, 



f( P ,C,VC) = fo(p,C) + -\VC\ 2 



(3) 



Here, /o is called the classical part of the free energy, 
C is the solute concentration, and e is a constant called 
the capillary coefficient. Coefficient e is assumed to be 
very small, so the second term is not-negligible only at 
the places of strong gradients of concentration, i.e. at 
interfaces. 

The chemical potential p, derived from the free energy 
function ([3]) with the use of assumption of incompress- 
ibility of mixture components (i.e. p is function of con- 
centration C and independent of pressure p) reads [3J 



p( P ,C) = /io(C) 



P dp 
p 2 dC 



-V • (pVC), (4) 



MC) = ^jP. (5) 



Here (1q stands for the classical part; the incompress- 
ibility assumption brings an explicit dependence of the 
chemical potential on pressure; and the last term is a 
non-classical contribution which stems from the Cahn- 
Hilliard addition in ([3]). 

The expression to be adopted for the classical part of 
the free energy function, /q , is written so as to reproduce 
the behaviour defined by the phase diagram depicted in 
figure [U In general, there are other different types of 
phase diagrams characterising the states of binary sys- 
tems [l7| . The diagram depicted in figure [1] is however 
one of the most popular. This diagram defines the mix- 
ture with the upper critical solution temperature: a sys- 
tem with the temperature below the critical point may be 



either homogeneous or heterogeneous depending on con- 
centration (the amount of a second component), while a 
supercritical system is always homogeneous. 

In the current work, we use the expression for the free 
energy function originally proposed by Landau [17j | to 
define a thermodynamic state of a system near its critical 
point, 



f (C) =a(C-C cr ) 2 + b(C-C cr ) 4 



(6) 



In this expression, C cr is the solute concentration in the 
critical point, and coefficients a and b are the phenomeno- 
logical parameters which define the choice of a particular 
binary mixture. It is necessary to note that the coeffi- 
cient a absorbs the factor (T—T cr ), i.e. (i) a tends to zero 
as the system approaches the critical point, and (ii) a is 
negative for undercritical conditions and positive for the 
temperatures above the critical value. Function (|B]) has 
two minima for negative values of a and one minimum for 
positive a. An equilibrium of a thermodynamic system 
corresponds to a minimum of the free energy function. 
Two minima characterising a system below the critical 
point are associated with two different phases, while in 
supercritical conditions a binary mixture is homogeneous 
and is characterised by the only minimum of the free en- 
ergy function. 

Using function ^ we may derive an expression for the 
classical part of the chemical potential, 



/xo(C) = 2a(C - C cr ) + 4b(C - C cr ) 3 



(7) 



If, firstly, the surface tension is not taken into account 
(the Maxwell theory is adopted) and, secondly, there 
is no loss or gain in volume upon liquids' mixing (sim- 
ple mixtures are considered), then the equilibrium con- 
centrations of the mixture components are defined by 
Coi,02 = C cr ± (— ~§b) 1 ^ 2 ■ 

We will use expression (JB]) as a simple model for a bi- 
nary system that, under different conditions, may define 
both miscible and immiscible solutions. Another popu- 
lar choice for the free energy function frequently used to 
model the behaviour of immiscible mixtures reads 



h = \{C-C m f{C-C a2 f. 



(8) 



It can be shown that ([5]) transforms into © if b = 1/2 
and Coi,o2 = C cr ± (—a) 1 / 2 . That is, formula ([BJ can be 
used to define the binary mixtures even under tempera- 
tures far from the critical one, where it can be considered 
as an interpolation of the experimental data; two new 
phenomenological coefficients a and b are determined so 
as to provide the best fit to the experimental data. 

Finally, we notice that a thermodynamic model de- 
fines the equilibrium states. A thermodynamic sys- 
tem can reach its equilibrium following a long equili- 
bration process (complete dissolution takes hours for the 
honey/water mixture). In the present work we aim to 
define the behaviour of a binary system during equilibra- 
tion to its thermodynamically stable state. 



4 



C. Hydrodynamic model 

The evolution of a binary mixture to its thermo- 
dynamic equilibrium is defined by the hydrodynamic 
model. The governing equations for hydrodynamic evo- 
lution of the Cahn-Hilliard fluid were derived in Ref. H. 
These equations include the laws of conservation of mass, 
species and momentum, 



| + v.(H = o, 



dC 

n \ — + (v ■ V)C ] = aV>, 



dv 
~dt 



(v ■ V)v = -Vp 



V • t v — eV • r £ 



(9) 



(10) 



(11) 



Here, the constant a is the coefficient of mobility. 

The fluids studied in this paper are assumed to be in- 
compressible, but due to large concentration gradients 
at the interface the continuity equation © needs to 
be taken in its full form. Such fluids are called quasi- 
incompressible [H,Q. The dependence of the fluid density 
p on the solute concentration C has the simplest form for 
the so-called simple mixtures, for which 



1-C 
Poi 



C_ 

P02 



(12) 



where p i and P02 are the solvent and solute densities. 
So for the mixture of water and glycerol, expression (|12[) 
is true with 1% accuracy; for the mixture of methanol 
and water, the accuracy of (fT2j) is 3.5% [l|. 

As the velocity field is non-solenoidal, the divergence 
term in the expression for the viscous stress tensor r n 
also needs to be taken into account, 



•„ = ?7 V (8> v + (V <g> t>) 



;(v-«)j 



(13) 



In this expression, the symbol <8> stands for the tensor 
product, the index T denotes the transposed, and / de- 
notes the unit tensor. The viscosity coefficient r], in- 
troduced by expression (flU)) , should, in general, depend 
on the solute concentration, and can be represented by 
a simple mass-weighted approximation or by an experi- 
mental formula. 

The tensor r e , defined as 



t £ = pVC (g> vc, 



(14) 



introduces the Korteweg interfacial stress. 

The right hand side of equation (|10l) takes into account 
the species transport due to diffusion. Here, the driving 
force for the diffusion flux is the gradient of chemical 
potential. The classical Fick's law which states a linear 
proportionality between a diffusive flux and concentra- 
tion gradient is only valid for weak solutions with small 



concentration gradients and cannot be used in our case 
One should also notice, that diffusive term of equa- 
tion (fTOj) includes the pressure and surface tension effects 
(stemming from expression for the chemical poten- 
tial). 

The full system of equations, defined in this section, 
is complemented with boundary and initial conditions. 
As initial conditions, the velocity v and concentration 
C are assumed to be known. As boundary conditions 
for the velocity field, the standard no-slip condition is 
used. For the field of solute concentration, two conditions 
need to be imposed at each boundary as ^ and (flU)) 
define the fourth order system of equations in terms of 
concentration field C. The first boundary condition is 
the zero-flux requirement, 



n ■ V/i = 0. 



(15) 



Here n is the vector normal to the boundary. The second 
boundary condition imposes the local thermodynamic 
equilibrium of the fluid-wall system. The full form of 
this condition is discussed in the works of [111 ] and [18| . 
For the simplest case of the contact angle of 90°, this 
condition can be written as 



n- VC = 0. 



(16) 



In such a form, this condition states that the interfacial 
boundary is orthogonal to the wall, i.e. molecules of the 
wall are neutral to components of a binary mixture. 

At the end of this section, let us enumerate the essen- 
tial features required for derivation of governing equa- 
tions ([9]), (fTUl) and (fTTj) . The first two requirements 
have been already mentioned, they are the use of the 
Cahn-Hilliard free energy function ([3]) and the assump- 
tion that the mixture components are incompressible liq- 
uids. The third requirement is to define the velocity, 
v, as the mass-averaged quantity over random veloci- 
ties of different molecules that constitute a fluid parti- 
cle. The volume-averaged definition for the fluid veloc- 
ity of a multiphase system is also possible and can even 
have some advantages as it would automatically produce 
a divergence- free velocity field for the mixing of two sim- 
ple incompressible liquids [l3, 19]. However, a traditional 
definition for the velocity in fluid mechanics is the mass- 
averaged quantity. The resultant governing equations 
based on this definition can be more easily and more nat- 
urally compared with the classical equations for a single 
phase medium. 



III. SEPARATION OF TIME-SCALES 



The quasi-compressibility of equations ©, (|10l) and 
(|11[) brings short-term processes. This significantly com- 
plicates the numerical solution of equations. The slow 
diffusion and convective evolution of a binary system 
would frequently be a primary practical interest, and 
even if we are not interested in quick processes, the nu- 
merical integration with a time step smaller than the 



5 



smallest typical time scale is required in order to ob- 
tain a convergent numerical solution. In this section, 
the time-scales that characterise the evolution of a mul- 
tiphase system are first identified and then, by using the 
multiple-scale method, the general governing equations 
are split into the systems which separately define the 
physical processes on different time-scales. 



The last term in equation (|2Tj) is the gravitational force 
(with 7 being a unit vector directed upwards). A non- 
dimensionalised coefficient of viscosity r\ (in formula (|25pl 
is obtained using the typical value rj* (e.g. the viscosity 
coefficient of a binary mixture in its critical point). 

The non-dimensional parameters entering the above 
equations are the Peclet number 



A. Non-dimensionalization of governing equations 

For further analysis we re-define the density and con- 
centration fields by shifting their reference points by the 
critical values, p cr and C cr , namely, 

(P - Per) -»• P, (C - C cr ) -> C. (17) 

We also non-dimensionalize the governing equations us- 
ing the following units of time r*, velocity V*, pressure 
and specific free energy: 

r* = V*=fx¥ 2 , p*=,o*^*, /*=/i*. (18) 

v # 

Here is the typical size, p* is the typical density (e.g. 
Per), and /i* is the unit of the chemical potential, for 
which we can use /x* = b. The resultant dimensionless 
equations then read 



^ + (».V)p = -(l + p)V-t>, 



<3u 



Pe 



(19) 



(20) 



(1 + P) ( gjj + ■ V)w 



-Vj> + — V • t„ - CaV • t 6 - Ga(l + p)7, (21) 



P= -(-M + Mo-Ca[AC + 0(l+p)(VC) 2 ]), (22) 



jo = AC 2 + PC 4 , M o = f^, 



1-0C" 



(23) 



(24) 



Pe = 



p*L* 



Otfl 



1/2 : 



the capillarity parameter 

Ca = 



P*L>1 : 



the Galileo number 



Ga = 



and the Reynolds number 



Re 



1/2 , 

v* 



(27) 



(28) 



(29) 



(30) 



The binary mixture is also characterised by the param- 
eter 



P* 
Poi 



p* 
P02 ' 



(31) 



where poi and po2 are the densities of the mixture com- 
ponents. The adopted expression for the free energy 
function also contains another non-dimensional param- 
eter that defines a distance from the critical point, 



P* 



(32) 



This parameter also defines whether the system is het- 
erogeneous (negative A) or homogeneous (positive A) in 
equilibrium. 

Finally, let us also write the expression for the total 
energy of the fluid enclosed by volume V, 



E= / (1 

Jv 



P) Y + /o + — (VC) 2 + Ga( 7 -r) dV. 

. ( 33 ) 

Here the first term corresponds to the fluid kinetic en- 
ergy, the second term is the classical part of free energy, 
the third term accounts for the energy accumulated by 
interfaces and the last term is the potential energy due 
to gravity. The total energy is measured in the units of 
p*p*L^. 



T r , = V ( V® v + (V® v) 1 - -(V • v)I ) , (25) 



B. Typical time scales 



Te = (l + P )VC ® vc. 



(26) 



To identify the different time scales characterising the 
physical system defined by the full equations written in 
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the previous section we consider an evolution of a one- 
dimensional (along the x-axis) small-amplitude pertur- 
bation on the background of a stationary homogeneous 
state. For simplicity, we also omit the gravity term for 
this analysis. 

Evolution of such a perturbation is defined by the fol- 
lowing linearised equations 



dp dv 
dt dx' 
dC 1 d 2 p 



1 

P = ~ 



dv 

at = 

2AC - Ca 



dt 
dp 
dx 



d 2 C 
dx 2 



Pe dx 2 ' 
4 d 2 v 
iRe^x 2 ' 

p = <f>C. 



(34) 
(35) 
(36) 
(37) 



Seeking a solution in the form of a plane wave, C ~ 
exp(ikx — iujt), where i is the imaginary unit, k is the 
wave-length and uj is the frequency, we obtain the follow- 
ing dispersion relation: 

^ = -iu{^ + ^e) + f 2 i2A + Can (38) 

Next, we split to into the real and imaginary parts, w = 
uj r + iuji and equate the coefficients in front of the like 
terms to obtain 

H - *?) - 



*W = -Ur (J + ^fc 2 ) ■ (40) 

The analysis of the second equation shows that this is 
satisfied by either 



LU r = 0, 



or 



1 (Pe 4, 2 



To simplify further derivations we will use that 
^ 2 ^^ 2 ) 1/2 «(5 + ^ 2 ) ; 



(41) 



(42) 



(43) 



or, equivalently, 

4> <C 1, A <C 1, Ca <C 1, Pe ~ 1. (44) 

That is, we assume a small density difference between 
components of binary mixtures, closeness to the critical 
point and the existence of an interface. 

Substitution of condition (|4ip into (j3"9")l produces two 
formulae for the decrements which define the monotonic 



decay of the considered ID perturbation, 

'Pe 4 



1 - • w7 lr 

^(2A + Cak 2 ) 



0JL2 



Pe , 4 12 
<l> 2 ~ 3i?e 



(45) 
(46) 



Next, we should note that equation (j3"9")l has no solu- 
tion if option (|4*2"j) to satisfy the imaginary part of the 
dispersion relation is chosen. 

Derived expressions (|4"5j) - (|46p allow us to conclude that 
the evolution of a multiphase system defined by equa- 
tions (|19p , (|20)l and (|2"Tj) is characterised by three differ- 
ent time-scales, 



Pe 



Td 



k 2 (2A + Cak 2 ) 7 



3Re 
Ik 2 ' 



Pe' 



(47) 



Here, Td and r c are the diffusion and convection time- 
scales, and r e is the fast time scale that we call the ex- 
pansion time scale. 

As it is clear now, the inequality (|43|) reflects the ra- 
tio between the quick and slow time-scales. For further 
analysis we will introduce the small parameter \ defined 
as 



X 



Te 
Td 



c/) 2 k 2 (2A + Cak 2 ) 
~Pe~ 2 ' 



(48) 



We are interested in dissolution dynamics, i.e. in evolu- 
tion of the system on the slow diffusive time-scale. 

The main idea of the multiple time-scale method is to 
represent the time derivative as 



d_ 

dt 



I d_ 

x 2 9t_ : 



X 



d_ 

dU 



(49) 



i.e. to explicitly introduce two times, the quick one, i_2, 
describing evolution of a system on the fast expansion 
time-scale and the slow one, t2, defining the diffusive and 
convective evolution. 

We will also expand all variables in the series of small 
parameter x^ 



X 2 V2 



4 

X V4 



P = Pa + X 2 P2 + X A Pi + 
P = XMi + X 3 P3 + 
P = X 2 P2 + X A Pi + 

C = X Ci+ x 3 c 3 + 
ri = l + 



(50) 
(51) 
(52) 
(53) 
(54) 
(55) 



and relate the non-dimensional parameters to different 
orders of x, 

<P = <PiX, A = A 2X 2 , Ca = Ca 2X \ (56) 

A, Pe = Pe , Ga = Ga 2 x 2 - (57) 
X 2 



Re = Re, 



We should note that there are two main conditions imply- 
ing the choices for the above- written ratios, namely, (i) to 
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save all essential physical effects and (ii) to have a closed 
system of governing equations. The non-dimensional 
parameters define the ratios between different terms in 
the governing equations. For derivations, we always as- 
sume that all effects in the original equations (fT9|) . ([20]) 
and (f2Tj) may be important and such values of the non- 
dimensional parameters are chosen which permit us to in- 
clude all corresponding terms in the final equations. Ra- 
tios (1561) imply that t ( i and r c are of the same order. This 
only means that the final equations will include both con- 
vective and diffusive terms. Following such an approach, 
we will obtain a comprehensive theoretical model that 
can be used for the analysis of a wide range of problems. 
For a particular problem, some terms in the model may 
be found either small or very large, which would, some- 
times, lead to further simplifications. But our aim is to 
obtain a general model. 



C. Separation of the processes occuring on 
different time scales 



First, we write down the different orders of equations 
(fT9|) . (pO]) and ([21]). The first orders of the mass ([19]) and 
species balances ([2"U]) are as follows 



° P2 =0, |^=0; 



at 



dt- 2 

i 



(58) 



9 Pa „ dC 3 i 

oT ± = -V-t; a , — 5- = — VV; 59 

OT-2 OT-2 Peo 

dC 5 ,dd ^„ dC 3 

Ot 2 Ot-2 



1 



VV (60) 



The first orders of the equation of momentum balance 
(EH) read 



dv 2 
~dt~2 



= -Vpo, (61) 



dv 4 



= -Vp 2 - G027, (62) 



dv 6 dv 2 dv A 
__ + _ + ( . 2 . V)t;2+p2 __ 



Re- 



-T n ,2 - Ca 2 r e ,2 - Ga 2 p2l- (63) 



Here, the corresponding orders of the viscous stress ten- 
sor and of the Korteweg tensor are 

2 

T v ,2 = V(g)t> 2 + (V® v 2 ) T - -(V -v 2 )I, (64) 
r e , 2 = VCiOVCi. (65) 

In these equations, all variables are assumed to be func- 
tions of both times, t_ 2 an d t 2 . Next, we will split out 
the processes on different time-scales using the averaging 
procedure briefly outlined in the next two paragraphs. 



Firstly, we assume that all fields can be split into 
slowly- and quickly-changing parts, 



V = u(t 2 ) + W(t-2,h), 

P = p(h) +p(t-2,h)- 



(66) 
(67) 



Here, u = ^- f^ 2 vdt and w = v — u. That is, u is 
the fluid velocity averaged over long time-scale (T 2 is a 
time interval much larger than the fast time scale) and w 
defines the quick-time-scale fluctuations of the fluid ve- 
locity. For scalar quantities, the barred symbol is used to 
denote the averaged parts and the tilded symbol denotes 
the fluctuating parts. As an example, we show the split- 
ting of the density field ([67]) , similar expressions have to 
be written for concentration, pressure and chemical po- 
tential. 

The equations for the long-term evolution will be ob- 
tained by averaging the equations. To accomplish aver- 
aging, the following general equalities are required to be 
used, 



V = 0, 



8 



dt- 



■(...) = o, 



(68) 



where V stands for any quantity. 

Averaging gives the following equations for the pro- 
cesses on diffusive and convective time-scales, 



Po = 0, fix = 0; (69) 
Vp 2 - Ga 2l = 0, (70) 



du 2 



■ + («2 • V)lt2 + {w 2 ■ \7)w 2 = -Vp4 

-V 2 u 2 - Ca 2 VCi <g> VCi - Ga 2 p 2 j, (71) 



V • u 2 = 0, (72) 
-^ + (u 2 .V)C 1 = —V 2 p 3 , (73) 

P2 = -?-(-£a + 2A 2 d + 4C? - Ca 2 V 2 d), (74) 

P2 = <f>iC 1 . (75) 

Here, we took into account that C% and P2 are indepen- 
dent of the quick time (|58[l so the over-bars were omitted 
for these variables. 

Subtracting the averaged parts from the full equations, 
we obtain the equations for the processes on the quick 
time-scale, 



dpi 
dt-: 



= -V • W 2 ; pA_=4>lCz, 



dC 3 1 



v 2 Ai 



ot_ 2 Pe 



dw 2 „ 1 , . . 

—_ = -Vp , Po = — (-Mi . 

OT_2 01 



(76) 

(77) 
(78) 



These equations describe a rapidly decaying process. Re- 
writing these equations in terms of one variable (e.g. 
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pressure), we find 



V 2 p 



V 2 ^ 



PO,in exp 



Peo, 



(79) 



where index 'in! denotes the initial value. If the aver- 
aged equations are characterised by the divergence-free 
fluid velocity (|72|) . for the quick processes, the quasi- 
comprcssibility effects are essential. But the divergence 
of the velocity field also rapidly decreases following a sim- 
ilar exponential decay, 



Omitting indexes and bars, the governing equations for 
the processes on the convective time scale can be finally 
written as follows 

(JUL I 

— + (u-V)u = -VII+— V 2 w-CV/z-Ga(/>C7, (81) 
at Re 



f + ( „.v) C = ijW 



V-it = 0, 



(82) 



(83) 



(V • w 2 ) = (V • w 2! , 



Oexp ( --s-t-2 



(80) 



Finally, we note that, firstly, the quick time-scale pro- 
cesses are damping and, secondly, there is no energy 
injection into fast processes as it could, for example, 
happen in the case of imposed high-frequency vibrations 
|21| . These two facts allow us to draw a conclusion 
that the short-term processes, even if existed at the ini- 
tial moment, are to be rapidly damped out and are not 
to affect the long-term evolution on the convective and 
diffusive time-scales. 

I 



2AC + 4C 3 - CaV 2 C. 



(84) 



Here variable II stands for the modified pressure, value 
of which can be determined using an incompressibility 
constraint. 

The derived equations must be supplemented with the 
following boundary conditions: 



v = 0, n ■ VC = 0, n • V/i = 0. 



(85) 



Finally, in order to show that no important effects have 
been lost during our derivation, let us write down the new 
expression for the total fluid energy: 



E 



u 

T 



+ /o + ?(VC) 2 + Ga0C( 7 ■ r) 



dV, f = AC 2 + C 4 . 



(86) 



The finally derived equations ([81]). (|82|) and (|83l) turn 
out to coincide with the model first used by [lCt 1 1 it] as a 
computational tool for tracking of the complex transfor- 
mations of an interfacial boundary between two immis- 
cible liquids. Expression © was used to define the free 
energy function. The computational model, adopted by 
[13 EL was earlier proposed by several researchers (for 
derivation and further references see Refs. [22| and HI). 
Here, we strictly showed that the equations of Jacqmin 
do, in fact, represent the Boussinesq approximation of the 
Cahn-Hilliard-Navier-Stokes equations. That is, these 
equations can be also used as the hydrodynamic model 
for slow dissolution processes in miscible systems. 



IV. CONCLUSIONS 

In the current paper, firstly, it was shown that the 
evolution of a multiphase system of two incompressible 
slowly miscible liquids is characterised by three typical 
time scales, namely, the convective and diffusion times 
and the typical expansion time (|47p . By using the mul- 
tiple scale method, we filtered out the short-term expan- 
sion process and, as a result, the equations for the slow 



I 

evolution of a binary system were derived. Even if the 
quick processes are not supported and should not affect 
the real physical system, the numerical integration of the 
full Cahn-Hilliard-Navier-Stokes equations would require 
time-resolution of all processes (the quick processes will 
be constantly excited due to discretization errors of a 
numerical scheme) and this would make calculations of 
long-term dissolution dynamics unfeasible. 

The obtained equations ([5T]). ([52} and (|53j! form a 
closed mathematical model for a general thermo- and 
hydro-dynamic evolution of a multiphase system. We 
should underline, that an applicability of the final model 
is to be considered in a broader extent compared to what 
is defined by assumptions (|56l) . and not, e.g., only for the 
binary systems in the vicinity of the critical point. For 
illustration, we may refer the classical Boussinesq equa- 
tions for thermal convection: despite the fact that such 
equations are derived in assumption of finite Rayleigh 
numbers, most interesting convection problems are con- 
sidered for large Rayleigh numbers (> 1000). 

We should also note that Refs. d, [ltj [HI and [H are 
all well known within the scientific community interested 
in hydrodynamics of multiphase flows. However, these 
papers contain different sets of governing equations mak- 
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ing other researchers to provide their additional justifi- 
cations for the actual models chosen, see e.g. Refs. 3 
and [H. Most of the numerical studies, based on the 
phase-field approach, use the Jacqmin's model but, in 
such papers, the evolution of immiscible systems is tar- 
geted. For miscible systems, the most popular approach 
is the set of impurity-like equations with an addition of 
the Korteweg stress tensor [4]-|6| . The main achievement 
of our work is an establishment of the relation between 
the full Cahn-Hilliard-Navier-Stokes equations of Lowen- 
grub and Truskinovsky and the divergence-free equations 
of Jacqmin. We showed that the Jacqmin's equations de- 
fine the thermo- and hydrodynamic evolution of both im- 
miscible and miscible systems. An important difference 
of Jacqmin's model from the impurity-like equations lies 
in the definition of diffusion flux through the gradient of 
the chemical potential rather than the concentration gra- 
dient. Such an amendment takes into account the surface 
tension effects in calculations of the dissolution rate. This 
means that the surface tension effects define not only the 
morphology of the interfacial boundary (the Korteweg 
stress in the Navier- Stokes equation) but also the diffu- 
sion rate. While the first effect is frequently taken into 
account, the second one is frequently missed out. 

We also would like to note that the derived equations 
dHU, ((52)) and ([55)1 are not ready yet for a complete anal- 
ysis of a particular configuration, as the values of the 
introduced non-dimensional parameters, e, a and b, are 
unknown. In fact, the value of a can be relatively eas- 
ily estimated from the phase diagram for a particular 
system. Determination of e and b would require a more 



lengthy investigation. Jacqmin published his equations 
about ten years ago. Since then, these equations were ap- 
plied to different problems |10l - [l4| , but they always were 
used for immiscible systems, for which a fluid behaviour 
in the limits of e — > and b — > oo is sought. Technically, 
the equations are successively solved for several values of 
parameters to reproduce the needed limiting behaviour of 
an immiscible interface, see e.g. Refs. [HJ El, El and [3 
For the case of miscible systems, the values of e and b are 
finite and should be obtained from a comparison of the 
numerical solution with the experimental data. In some 
recent experiments, the estimations and measurements 
of the surface tension coefficients for several particular 
miscible binary systems became available [24T - [2q | . These 
experiments also contain detailed description of the dis- 
solution dynamics. Such information (time evolution of 
the surface tension coefficient and of the dissolution rate) 
should be sufficient for obtaining the missing values of e 
and b. This is the aim of the author's current research 
work. 
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